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Hopping transport in a one-dimensional system is studied numerically. A fast algorithm is devised 
to find the lowest-resistance path at arbitrary electric field. Probability distribution functions of 
individual resistances on the path and the net resistance are calculated and fitted to compact 
analytic formulas. Qualitative differences between statistics of resistance fluctuations in Ohmic and 
non-Ohmic regimes are elucidated. The results are compared with prior theoretical and experimental 
work on the subject. 
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I. INTRODUCTION 

It is well known that low-temperature transport in dis- 
ordered one-dimensional (ID) structures is distinguished 
by large mesoscopic fluctuations. Such fluctuations have 
been measured^ 2345 gyg^ samples of considerable 
length. They arise from the interplay of localization and 
rigid geometrical constraints on possible current paths. 
The total resistance tends to be dominated by a few 
strong obstacles — "breaks" — - which occ ur at random 
due to disorder in the sample.^* ^ * ^ * ^ * ^° * ^Q This unusual 
behavior can be contrasted with a more familiar case of 
dimensions d > 1. There the current can go around the 
breaks, so that the mesoscopic fluctuations of transport 
properties are usually small and self-averaging. 

In this paper we consider ID systems that are not 
too short, so that the c oheren t tunneling of electrons 
through their entire lengtbP^'^is extremely improbable. 
Instead, electrons traverse each sample via a sequence 
of many incoherent tunneling acts — the variable-range 
hopping^ ( VRH) . By studying the VRH transporili^i one 
aims to extract information about the nature of electron 
localization and disorder in the system. However, this 
task is far from trivial. Although the basic physics of 
the ID VRH problem is quite well understood, experi- 
mental studies of VRH are typically done in a narrow 
parameter range where usual theoretical approximations 
are still rather crude. Below we demonstrate that large 
corrections appear when the transport properties of a 
standard VRH model are calculated numerically, which 
means, with fewer approximations. 

To deal with large mesoscopic fluctuations we follow 
prior work a nd comput e both the probability distribu- 
tion function^iSmilllill (PDF) and suitable averages of 
the transport observables. For example, we study the 
ensemble-averaged conductance (G), which can be mea- 
sured experimentally by having a large number of ID 
wires connected in parallel.'^ 

Our primary purpose is to investigate non-Ohmic ef- 
fects, e.g., the dependence of function {G){F,T) on the 
electric force F — —eE. This regime has been studied 
much less compared to the Ohmic one. However, re- 
cently an analytical theory of non-Ohmic ID VRH has 
been proposed in a work of one of us. 1^ Here we approach 



the same problem numerically. We have developed an ef- 
ficient computer algorithm, which is able to find the VRH 
conductance of a given sample at arbitrary electric field. 
By choosing a low F the Ohmic conductance G(0, T) can 
also be calculated. 

Since the Ohmic case has been more widely studied, 
it deserves a brief discussion first. In dimensions d > 1 
the Ohmic conductance is known to follow the stretched 
exponential temperature dependence: 

G(0,r) = Goexp[-(A/r)^] , (1) 

where A is some energy scale, and Go depends on T 
at most algebraically. The exponent 7 = l/(d+ 1) at 
d > 1 signifies the Mott law. The Efros-Shklovskii law 
corresponds to 7 = 1/2. It applies when the long-range 
Coulomb interactions are important .^Sl 

In ID, the Mott and Efros-Shklovskii exponents coin- 
cide. This is because in ID the 1/r Coulomb potential 
is only marginally long-range to begin with, and then 
typically also screened by a nearby metallic gate. The 
importance of the remaining interactions is determined 
by the dimcnsionless parameter 

e = 1 + (e'g/C) , (2) 

which has the physical meaning of the dielectric constant. 
Here C is the capacitance to the gate per unit length of 
the wire and g is the average density of states. (Note 
that e is related to the Luttinger-liquid parameteiS^ of 
a disorder-free ID system.) In this paper we study the 
case of weak interactions, e ~ 1, where, naively, the Mott 
law may seem to be a reasonable starting point. 

Actually, the ID Mott law is modified by the afore- 
mentioned mesoscopic fluctuations. Le^ and Raikh and 
Ruzin^^ showed analytically that at low temperatures the 
energy scale A in Eq. ([T|) is not a constant but a logarith- 
mic function of T. More importantly, A(T) is determined 
not only by intrinsic properties of the system but also by 
its size. As T increases, a narrow range of temperature 
appears where another dependence, A(T) cx 1/T is real- 
ized. Hence, instead of the 7=1/2 Mott law we effec- 
tively have a simple activation,!^ 7 = 1. Su ch behavior 
has been conflrmed by numerical simulationsJ ^ * ^° * ^^ * ^^ * ^° l 
so it is considered well established. 
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Nevertheless, to establish a reference point for our 
study of non-Ohmic VRH we examined the Ohmic con- 
ductance carefully by our method. Remarkably, we found 
that it is essential to introduce often discarded "sublead- 
ing" terms in the analytical expressions. If this is not 
done, analytical and numerical results for G can differ by 
orders of magnitude. 

In the non-Ohmic regime, w hich is our main subject of 
interest, it has been customarjl^lllllllllll to characterize 
the field-dependence of the conductivity by means of the 
length parameter Lc'. 



{G)iF,T) ^ (G)iO,T)exp{\F\L,/T) 



(3) 



In experiment, this law typically describes the first 
decade of the conductivity rise . Ther eafter, deviations 
tend to occur. Indeed, in theorjEHHEZI jg expected to 
be not a constant but a function of F and T. We will 
show that in ID Lc may also depend on the averaging 
procedure utilized to obtain (G). 

At large enough F, Eq. ^ eventually becomes a poor 
approximation. Theoretically, it should cross over tcP 



G 



2LRo 
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where Tq is defined by 



To = 1 / (ga) 



(4) 



(5) 



a is the electron localization length, and Rq is specified 
in Sec. Ill (At such fields mesoscopic conductance fiuc- 



tuations are small, and so we denote (G) simply by G.) 
Our numerical results are consistent with Eq. Q. Note 
that it can be viewed as the ID Mott law with the ef- 
fective temperature^-'^ Tcff ~ Fa replacing the ambient 
temperature T. 

Finally, we examine the PDFs of the mesoscopic con- 
ductance fluctuations. Such functions can also be stud- 
ied experi ment ally, albeit it requires a substantial time 
and effort .122^1] We demonstrate that the PDFs are qual- 
itatively different in the Ohmic and non-Ohmic regimes. 
Both have asymmetric long tails. However, the Ohmic 
PDF is skewed towards the low conductances, while the 
non-Ohmic one towards the high conductances. We ex- 
plain these differences and show how they evolve as a 
function of the applied field F. 

The paper is organized as follows. In Sec.|ll]we present 
the summary of our results. In Sec. |III| we define the 
model and introduce our fast algorithm for computing 
the resistance at a given current. In Sec. IV we obtain 
analytical fitting formulas for the PDF of individual hops 
in the Ohmic and non-Ohmic regimes. We also describe 
approximate but much faster "PDF-algorithm" for com- 
puting the net resistances. Section |V] discusses the dif- 
ferences of two averaging procedures: at given current / 
and at given electric field F. Finally, Sec. |VI| contains 
discussion and comparison with experiments. 
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FIG. 1: Ensemble-averaged Ohmic conductance (G) as a func- 
tion of temperature (the curves with fluctuations): (a) Rel- 
atively high T. The dashed line is the best fit to the simple 
exponential law, 7 = 1 and A = 0.62ro in Eq. ([l]). (b) A range 
of low r. The dashed line is a fit to the ID Mott law, 7 = 1/2 
and A = 8.4To in Eq. The upper curve is Eq. The 
dots are the Ohmic limit of the upper four traces in Fig. [3] 



II. MAIN RESULTS 

In this section we provide a short overview of our 
principal results for experimentally measurable transport 
properties. 

Figure [T] shows the dependence of the average Ohmic 
conductivity (G(0,T)) on temperature in an ensemble 
of samples of length L = 250a. To test the expected 
crossover behavior, we fit the low T data points using 
Eq. ^ with 7 — 1/2, corresponding to the ID Mott 
law. We fit higher T using 7 = 1, representing activated 
transport. In the Mott regime we find A = 8ATq. For 
the activated regime we get A = 0.62To. Note the large 
difference between these values. As far as A is concerned, 
our numerical results are in a good agreement with the 
analytical theory of Raikh and RuzinSl (RR) . In the high- 
T regime it predicts A = To/2. Their low-T formula 
reads 



G = i?Q ^ exp 



To 
T 



(6) 



where 1/ is defined as the solution of the transcendental 
equation 



2T ( L 
V — in \Jv — 
To \ a 

Therefore, RR result for Mott's A is 



A(r) ^ 2To In ( y/iy - 
a 



(7) 



(8) 



Strictly speaking, it is not a constant but a slow function 
of T. In the range of T where the fit to the Mott law 
was done, it is indeed close to 8.4To. The large difference 
between the values of A in the Mott and the activated 
regime is due to the "large" logarithm ha{y^ L/a). 

When the RR formula is plotted alongside our numer- 
ical results, it is seen to exhibit a very similar functional 
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FIG. 2: The PDF of the logarithm of the total resistance 
R in the Ohmic limit. The simulations are done for system 
size L — 10^, localization length a — 4, and um = 12.247. 
The smooth curve on the right is obtained using the PDF 
algorithm; the markers correspond to the shortest-path sim- 
ulation. The leftmost curve is Eq. (l9|. 




FIG. 3: Conductance as a function of a scaled electric field 
Fa /To (five solid lines on the left). The simulations are done 
for system size L — 10^ a nd localization length a = 4. The 
values of um = \/2To/T are indicated next to each curve. 
The fits to Eq. ([3| used to extract Lc are shown by the dotted 
lines. The rightmost curve is Eq. (|4|. 



behavior yet a large difference in the absolute value, see 
Fig. [ijb). Despite the fact that we study exactly the 
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same model, see details in Sec. Ill RR's predictions dif- 
fer from our results by two orders of magnitude. We 
attribute this discrepancy to the "subleading" terms not 
included in the asymptotic theory of RR. 

Next, we present the PDF Pu{U) of the logarithm of 
the total resistance U = ln(i?/i?o) in the Ohmic limit. 
Fig. |2] for the same set of wires at temperature T ~ 
ro/75. Our curves are plotted side-by-side with RR's 
formula 

P[/(C/) = ^^exp[-^A;JC/-exp(-^<5^7)] , (9) 
SU = U-{V^To/T). (10) 

Again, we see that while the shapes of the curves are 
practically identical, RR's distribution is centered around 
a lower value of U. This is consistent with the difference 
of the G{T) curves described above: ignoring the "sub- 
leading" terms results in a decreased resistance. 

Let us now turn to the non-Ohmic regime. Figure [3] 
illustrates the dependence of the ensemble- averaged con- 
ductance as a function of the applied electric field at five 
different fixed T. At low fields the conductance strongly 
depends on T, as the curves originate at points on the 
vertical axis which differ by many orders of magnitude. 
[Four of these points are also shown as dots in Fig. [ijb) .] 
All the traces grow monotonically with F. Equation ^ 
gives an adequate fit (dotted lines) in a range of low 
fields. The corresponding Lc are presented in Fig. |4] 
We plot them as a function of both the temperature and 
the "Mott value" , , defined as 



UM = (2To/r)i/2. 



(11) 



We see that Lc ~ I.Oumo, which is the average hop 
length. This implies that the average conductance is 
dominated by rare samples that do not contain large 




FIG. 4: Characteristic length Lc [Eq. ([3|] that determines the 
non-Ohmic behavior as a function of temperature (dots). For 
comparison, the dashed curve represents the relation Lc/a = 
1.9itAf, which corresponds to a typical hop length. 



breaks, so that the total voltage is distributed roughly 
equally among all the hops. In contrast, we know that 
the average resistance is determined by typical samples 
where the breaks are present; the entire voltage is ap- 
plied to the single most resistive hop, and the size of the 
non-Ohmic effect is much larger, see Fig. |5] We discuss 
the difference between average conductance and average 
resistance in more detail in Sec. IVII 

At large F the rise of the conductance becomes less 
rapid than exponential and the curves in Fig. [3] tend to 
converge to a common T-independent envelope of Eq. (|4| , 
confirming the analytical predictions of Fogler and Kel- 
ley.'^ At such high electric fields F, high-resistance 
breaks are eliminated not only from rare samples but 
from typical ones. This can be deduced from the fact 
that averaging of the conductance G approaches the re- 
sult of averaging of the resistance R (followed by taking 
the inverse). As evident from Fig. [s] the two curves in- 
deed approach each other with increasing field. A de- 
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FIG. 5: Dependence of the conductance on the scaled electric 
field averaged in two different ways. The upper line is the av- 
erage conductance, the lower one is the inverse of the average 
resistance. Simulation parameters are the same as in Fig. [2] 



tailed analysis of this crossover in terms of the PDFs is 
given in Sec. |IV| 

This concludes the summary of our main results. In 
the next section we define the model and the method of 
calculation by which they have been obtained. 



III. MODEL 

A. VRH resistor network 

We model our samples as a network of resistors, as is 
customary in the VRH theory.^^^l To derive the parame- 
ters of this network we proceed as follows. The phonon- 
assisted transfer of electrons from one localized state (LS) 
i to another j is characterized by the transition rate 



ro/.(i-/,)x 



r iV(Ae) 



if Ae > , 



\ iV(|A£|) + 1 , otherwise , 



(12) 



where fi is the occupation factor of ith LS, N{e) is the 
Bose-Einstein distribution, and Ae is the energy differ- 
ence in the hop: 



Ae = Ei 



(13) 



Here e^ and e^ are the energy of i th LS with and without 
the applied field, respectively, and is its electrostatic 
potential shift. In a realistic model, the rate prefactor 
Fq should have some algebraic dependence on Ae, which 
counteracts the divergence of A^(|Ae|) at Ae 0. How- 
ever, such Ae are virtually never important in the VRH 
transport. For simplicity, we treat Fq as a constant. 
The net current between the LS i and j is given by 



(14) 



In order to compute lij, one needs to know the occu- 
pations factors of all LS. They can be found from the 



conditions of current conservation (the so-called Master 
equation) , 



(15) 



supplemented by suitable boundary conditions at the 
source and drain electrodes. Unfortunately, these equa- 
tions are nonlinear and involve an exponentially large 
spread of the values of fi. This makes the solution dif- 
ficult to obtain. It can be done nu merically, using some 
clever iterative techniques .I^SlSSISll However, the rate of 
convergence is slow. We proceed in a different direction, 
which enables us to map the problem to a resistor net- 
work even in the non-Ohmic regime. As a result, we can 
achieve practically the same speed of simulations in the 
non-Ohmic regime as in the Ohmic one. 

We start by defining the chemical and the electrochem- 
ical potentials as follows: 



fi^ = nn(/, ^ - 1) , 'qi = ^li- e<^i 



(16) 



The "voltage drop" of every (i, j) link is given by the 
difference of electrochemical potentials 5rj — rfi — rjj . In 
turn, the link resistance is defined by 

R,, = S^/I,j . (17) 

Substituting this into Eq. (14 1, one obtain^23l 



/=^sinhf-^ 
eRo \ 2T 



exp 



'2Xi 



X exp 



2T 



2T 



2T 



(18) 



where Xij is the distance between the LS i and j, and 
Ro = T/{e^To). 

Let us introduce logarithmic variables 



1 ^ij 

= In — = ui 



In 



Si] 

Y 



uj = In 



T 
^Rol 



(19) 



It is easy to see then that if the voltage drop is smaller 
than T (Ohmic case), the expression for u.^ reduces to 
the well-known form-'^^ 



2xi 



v\ 



2T 



2T 



2T 



(20) 



Here either 77^ or -qj can be used for 77. 

To complete the system of equations, we need a for- 
mula for the electrostatic potential <i>i [Eq. ([l3])]. It is 
determined by charges on the source and drain leads, 
and the perturbation of the electron density inside the 
wire (given by the occupation factors fi). The relative 
importance of these contributions depends on the exact 
geometry of the device. We consider a typical situation 
where there is a metallic gate positioned parallel to the 
wire, with C again denoting the capacitance to the gate 
per unit length of the wire. We further assume that the 
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capacitive coupling to the leads is much smaller and can 
be neglected. In this case, we find 



en{x) 



(21) 



where n{x) is the deviation of the local density from equi- 
librium. Neglecting fluctuations in the local density of 
states and any correlation effects, we can directly relate 
n{x) to the local chemical potential, n(x) = gfi(x), which 
implies 



-e<I>, 



2 

C 



e 



(22) 



where e is given by Eq. In comparison, m previous 
literature it was common to approximate <i>i simply by 
—Fxi, i.e., to assume that the electric field in the system 
is uniform. Although this may be reasonable for a sample 
of dimension d > 1 with bulk leads, it is inappropriate 
for the specified ID geometry where the electric field is 
heavily concentrated at the breaks. 



Substituting Eq. (22 I into Eq. (18 1, we obtain 



2T 



- sinh 



(577 \ 



X exp 



X exp 



2T 



2T 



+ (r;,-ry,-)(6-l)/e| 
2T 



(23) 



In this equation, all self-consistent field effects are conve- 
niently expressed in terms of the effective dielectric con- 
stant e. Actually, in this paper we focus on the case of 
weak electron interaction, so that henceforth e will be 
replaced by unity. Effect of finite-strength interactions, 
e 7^ 1, will be considered in a separate publication. 

To implement the resistor network we proceed as fol- 
lows. We choose the coordinates of the LS, < Xi < L, 
to be the sites of a chain with unit nearest-neighbor spac- 
ing. Their energies e° are selected randomly. We draw 
these energies from the Poisson distribution Pe{z) — 
gexp{—g\z\) and generate two of them — one above zero 
and the other below — at each internal lattice point. 
For high currents (small u/) we sometimes generate ad- 
ditional energies at the same lattice point, using the same 
procedure. 

The leftmost lattice point is the source electrode. It 
has only one LS at the coordinates (0,0), whereas the 
right end of the sample has many sites at the same ap- 
position, equally spaced along the energy axis, see Fig. |6] 
This is done in order to simulate the behavior of a metal- 
lic drain electrode where there are all energies present. 



B. Shortest- path algorithm 

At this point, we make a crucial approximation, which 
is, however, conventional in the VRH theory.^'' We will 
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FIG. 6: An example of the optimal path in a modestly non- 
Ohmic regime, uj — 25. The dots represent localized states. 



suppose that there exists a certain path through the net- 
work — the optimal path — whose conductance is much 
higher than any other linear path. We can assume then 
that all the current flows along the optimal path without 
branching. As we show below, this allows us to devise a 
fast algorithm for flnding such a path and therefore the 
net resistance of the sample. 

1, we can express the voltage 



Using Eq. ([18| with e 
drop Sr] in terms of 77^, the bare site energies and e^, 
and u/. To this end we define axillary variables t and q: 



i=(e"-77.)/r. 



-3 
2x 



Q = 



2T 



2T 



ui 



(24) 
(25) 



Only q < Q are physically allowed, which means that 
there is a certain maximum current that can fiow through 
the given link. If so, the voltage drop in question is 
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-ln(l-e«) , 
In (1 + e«~*) , 



if > 1 - 
otherwise. 



(26) 
(27) 



One can show that this cumbersome expression is reduced 
to the familiar Eq. (20) in the low-current limit, u/ 00. 
Indeed, in this case, Eq. (26 1 applies for i < 0, while 



Eq. (|27| for t > 0. Note that Eq. ([201 is independent of 
M/, as is appropriate in the Ohmic regime. 

We can use the above equations to find the optimal 
path through the sample. This is the path which would 
require the lowest voltage (difference in the electrochem- 
ical potential between the ends of the sample) for a given 
current. To do so we use the well-known Dijkstra algo- 
rithm^^ to calculate the minimum "cost" of getting from 
the source to the drain. Here the cost is the total volt- 
age V . Similarly, the cost Ci of getting to site i on the 
optimal path is 



(28) 



The algorithms starts by assigning zero cost to the source 
(0, 0) and infinite cost to all other sites. Thereafter the 
spanning tree of the lowest-cost sites is grown iteratively. 
Initially, the tree consists of only the source site. At each 
iteration, a site of the lowest cost among those that are 
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still outside the tree is added to the tree. The costs of 
sites j outside the tree are relaxed (updated) according 
to the rule 



(n+l) . / (ri) 

c , = mm [^c - 



St], 



(29) 



(n) 

Here c) is the cost of site i at nth iteration. The cost 



increment Srj in Eq. (29 1 is computed using Eqs. (26) and 
( 27 1 . The process terminates when any of the LS located 
on the drain electrode are reached. In Fig. [6] one can see 
an example of an optimal path found by our algorithm 
in a modestly non-Ohmic regime. 

In the Ohmic VRH problem, the Dijkstra algorithm 
has been used in Ref. [TS] In that regime each link has a 
fixed cost. Here we are using the Dijkstra algorithm in 
an unconventional situation where the cost Srj = Sr]{ci) 
of a given link is not a constant but a nonlinear function 
of the cost of the earlier sites in the tree. A potentially 
troublesome point is that in the course of iterations we 
retain only the lowest cost so far. We effectively assume 
that for any i and j 



mm c; 



(57/(min Ci) . (30) 



Let us show that this equation is satisfied, which implies 
that our algorithm works correctly at arbitrary current. 
First of all, by our earlier assumption the current does 
not branch, and so the current through any link of the 
optimal path must be exactly /. Second, a sufficient 
condition for vahdity of Eq. (30 1 is dcj/dci > 0. That 
is, increasing Ci by taking a less optimal path to the ith 
site would not help to decrease cj . In view of Eq. ( 28 1 , 
the last condition can be written as 

d 



5r]<l. 



(31) 



We need to examine the two possible cases represented 
by Eqs. (26) and (27). In the former, we get 



d 



Srj 



% ' 2(1 -e9) 
In the latter, we obtain 



[sgn(£, - ?7^) + 1] < < 1 . (32) 



drji ^'^ 2(1 + e9-*) 



— [l-sgn(e,;-77i)] < 1 



(33) 



In both cases inequality (31) is satisfied, which means 



that our algorithm does find the optimal path. 

In the course of simulations, the resistance of every 
link on this path as well as their total sum are saved for 
further analysis. Repeating the process over many dis- 
order realizations, we obtain the PDFs and the averages 
of desired transport properties, discussed in more detail 
below. 



IV. DISTRIBUTION FUNCTIONS 

In this section we review analytical predictions regard- 
ing the functional form of the PDF of link resistances 
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FIG. 7: Geometry of (a) Ohmic (b) Non-Ohmic break in the 
energy-position space. The dots represent localized states. 



and compare them with the simulation results. In both 
Ohmic and non-Ohmic cases we are able to make the 
two to agree by introducing a few refinements in the an- 
alytical formulas and by adjusting numerical coefficients 
therein. 



A. Ohmic case 

We start by discussing the Ohmic case: u/ oo. Ac- 
cording to previous theoretical studies, notably Refs. 121 
[TUl and [m the logarithm of the average resistance of a 
link is on the order of the Mott value um- Links with 
u S> Um are exponentially rare; however, they act as 
bottlenecks and the total resistance depends on them. In 
order for such high-resistance links to exist, the optimal 
path has to encounter regions in the energy-position (x-e) 
space that are empty of LS. Using the method of optimal 
fiuctuation, RRpi] showed that the leading asymptotic 
behavior of the PDF of the breaks has the form 



du 



exp[— gA(u)] 



(34) 



where A{u) is the smallest possible area of a break with 
given u in the x-e space. Equation (34 1 is due to the Pois- 
son distribution of the LS in the x-e space. The shape 
that attains the minimal area depends on whether the 
break is Ohmic {u < uj) or non-Ohmic {u — uj 3> 1). For 
the former case, RR showed that the break is diamond- 
shaped with the width ua/2 and the height 2uT, see 
Fig.|7|a). This entails the quadratic dependence 



gA{u) = {u/uAif 



(35) 



Later, taking into account shape fluctuations of the break 
along its perimeter, Ruzin'^^ proposed a refined formula 

P{u) = Coexp(2BM/uM) x gA' {u) eyi-p[-gA{u)]. (36) 

While Co is determined essentially by the normalization 
of P, analytical calculation of the coefficient B is chal- 
lenging. Ruzin gave a rough estimate B « \/2/3 « 0.5. 
In this study, we calculate B numerically. Indeed, from 
the example of the optimal path shown in Fig. |6j it is 
clear that the voids around the long hops hardly ever 
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look like "diamonds" (or "hexagons", see below). This 
means that even though the RR theory provides the basis 
for understanding the behavior of P{u), numerical simu- 
lations are critical in order to calculate it accurately. 

At each L the functional form of the P{u) is expected 
to depend only on the dimensional ratio u/um- By run- 
ning simulations at different combinations of a, g, and T, 
we convinced ourselves that this is indeed correct, for the 
exception of very small u where lattice discreteness starts 
to matter. Fortunately, such u are irrelevant for the 
macroscopic transport properties as they do not deter- 
mine the resistance. Thereafter we fixed a = A, g — 1/3, 
and T — 0.01, which yields the characteristic tempera- 
ture To = 3/4 and the Mott parameter um — 12.247, 
cf . Eq. (111. To ensure we are in the Ohmic regime 
ui = 200 ^ Um was used. 

For each L in the set L = 100,200,400,500, and 1000 
we generated many realizations of ID wires, respectively, 
20000, 10000, 5000, 4000, and 2000. Anticipating the 
finite-size effects, these numbers were chosen in order to 
have the same total number 2 x 10^ of LS at each L. We 
found optimal paths through the samples and created 
the PDFs of the link resistances. We fitted such PDFs to 
Eq. p6| using _B as a single adjustable parameter. The 
quality of the fits was rather good, see an example in 
Fig. [Sj Furthermore, even though Eq. (36 1 is meant to 
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apply at u S> um , it fits our numerical results for u < um 
as well. 

Interestingly, we found that B slowly but systemati- 
cally increases with L. When plotted as a function of 
l/L, it was seen to vary linearly, tending to a constant 
for large L. We believe that the reason for this finite-size 
effect is the following: due to the source electrode being 
at zero energy, the resistance of the first link is typically 
lower than average. In shorter samples, where the total 
number of hops through the sample is about ten or 
so [see Eq. (41 1 below], it impacts the PDF. As the sam- 
ples get longer, increases and this first hop does not 
influence the overall PDF any more. To get the value of 
coefficient B in the thermodynamic limit, we used linear 
extrapolation to L = cxd. Our final estimate is 



5 = 0.92 ±0.02, 



(37) 



approximately twice larger than that of Ref . 1361 

Two characteristic measures of the width of the distri- 
bution are its mode and its average. For P{u) they are 
given by, respectively. 



Um 



(1.30 ±0.02) MM, 



(38) 



(u) = j uP{u)du = (1.39 ± 0.02) um ■ (39) 



As expected, both are the order of the Mott parameter 
Um- One more important quantity is the average num- 
ber Nu of links on the path. It determines the relation 





FIG. 8: Numerical results for P{u) in the Ohmic regime 
shown on (a) linear and (b) logarithmic scale. The simulation 
parameters are the same as in Fig. [2] e.g., um ~ 12.247 (thin 
line). Th e small fluctuations are of statistical origin. Equa- 
tion (361 with B = 0.9 is represented by the smooth thick 
line. 



between P{u) and the probability density of breaks per 
unit length of the wire p{u) : 



p{u) 



L 



P{u) 



(40) 



Since the width of each link is not smaller than {a/2)u, 
cf. Eq. (lis]), Nu can be estimated from below as 
{2L/a)/{u)Ki lAL/auM ■ According to our simulations, 
the actual is approximately twice larger: 



Nu = (3.04 ± 0.07) 



Mm a 



(41) 



Besides RR— and Ruzin^^^, the calculation of P{u) was 
previously attempted by Ladieu and Bouchaud!-^ They 
reported Mmax and (m) that differ from our Eqs. (38 1 and 



(39 1 by 30-40% In fact, we were unable to verify that 



statement because the main equation of Ref. |37] has no 
solution. As written, that equation does not conserve 
probabilty. Consequently, we believe that our results 
constitute the first reliable calculation of function P{u). 



B. Non-Ohmic case 

Let us now discuss the breaks in the non-Ohmic 
regime. Unlike the diamonds of the Ohmic case, the non- 
Ohmic breaks are hexagonal, see Ref. 18 and Fig. |7] with 
area 



gA{u) 



w"^ + w(3 



w — uj + ln(l 



-/3 



'■M 



(42) 
(43) 



The width of the break in the real space is wa/2. Note 
that at u < uj we have wa/2 ~ Ma/2, which is the width 
of the Ohmic break. The combination f3T, which is equal 
to the electrochemical potential drop across the break, 
gives the the height of the middle part of the break in 
the x-e space. 
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In order to account for the possible perimeter correc- 
tions to P{u), we consider the following trial form: 



0.15 



P{u) — Co exp 

X gA'{u)e^^[-gA{u)] 



2B^ + C(^ 
um \ um 



(44) 



Here the contribution of the top and bottom parts of the 
perimeter is modeled after Eq. ( 36 1 . It is proportional 
to the length of such parts ~ w and the coefficient B. 
The contribution of the side walls of the break, of length 
/3T, is written differently. Indeed, Ruzin's argument'^^ 
suggests that they give no contribution at all. In fact, 
we found it necessary to include a correction albeit with 
a smaller exponent D — 0.5. We have no other justifi- 
cation for this exponent except that it provides a good 
fit to the numerical P{u), see below. The explicit for- 
mula for P{u) can be derived from Eqs. (42|-(44l by the 
straightforward differentiation with respect to u. How- 
ever, it is cumbersome and we do not write it here. Equa- 
tion (44) applies for u — u/ ^ 1 and ui ^ um- It refines 



the corresponding expression for P{u) in Ref. IF" where 
the first (subleading) exponential term was not included. 
The Ohmic and non-Ohmic formulas, Eqs. (36 1 and (44 1, 
match at u 



Ul 



1. 



Equation ( 44 ) predicts that P{u) decays as a Gaussian 



dX Um < u < uj and as an exponential of the exponen- 
tial at u > uj. In between, it exhibits a narrow peak 
of width Su ^ lii{u\j /uj) near the non-Ohmic threshold 
u — uj. For parameters chosen in Fig. |9] this peak is 
so pronounced that it already dwarfs the "Ohmic" maxi- 
mum at u — Mmax- The reason for its appearance is sim- 
ilar to that discussed in a three-dimensional case.'^^ This 
narrow peak is due so-called "soft" links that used to have 
resistances u > m/ in the Ohmic regime. Such links are 
similar to forward-biased diodes: their conductance in- 
creases exponentially with the electrochemical potential 
drop Srj. When a finite current is made to fiow across 
the wire, such links self-generate Srj large enough to push 
their resistance back to an immediate vicinity of the non- 
Ohmic threshold u uj. 

The soft links are realized when the energies at their 
endpoints satisfy a certain inequality, which can be de- 
rived from Eq. ([Ts]) or looked up in Table I of Ref. [231 
Therefore, not all links are soft. There also also "hard" 
links, which are similar to reverse-biased diodes, whose 
resistance does not change much with Srj. These links 
are never included in the optimal path because they are 
simply not able to support the necessary current J. The 
peculiar shape of P{u) that follows from these arguments 
is nicely confirmed by simulations, which we now briefiy 
describe. 

The simulation procedure in the non-Ohmic regime 
is practically identical to the Ohmic case with one ex- 
ception: we have to put more than two energy sites at 
each lattice point x. The reason for this is that for high 
currents (and, therefore, high voltages), as the electron 
moves through the sample, it hops onto LS with lower 





FIG. 9: Numerical results for P{u) for finite current, uj = 
20, shown on (a) linear and (b) logarithmic scale (thin line). 
The small fiuctuations are of statistical origin. The fitting 
formula (|44| with um = 12.247, B = 0.9, C = 0.75, and 
D = 0.5 is represented by the thick line. 



energies, see Fig. [6j The addition of extra LS is done 
to ensure that there are LS for the electron to hop onto, 
otherwise the path would not be found. We also have 
to increase the range of the energies on the electrode for 
exactly the same reason. The simulation was conducted 
at uj = 35,30,25, and 20. Two values of um are used: 
12.247 (same as above) and 20 (obtained by adjusting 
the temperature but keeping g — l/S the same). The fit 
of the numerical P{u) to Eq. (|44| for um = 12.247 can 
be seen in Fig. |9] and it is quite good at all but very small 
u (which are irrelevant, see the note above). 



C. Distribution of the net resistance 

Besides studying the distribution of individual hops, 
we also investigated the statistics of the net resistance R. 
In Fig.[TO]we present the a sequence of four PDF's of C/ = 
ln(i?/i?o) obtained from our shortest-path simulations. 
From one curve to the next the current increases by the 
same factor of exp(5). A qualitative difference from the 
PDF for the Ohmic case (Fig.|2| is immediately apparent. 
The Ohmic PDF is skewed to the right, towards the large 
resistances. In contrast, the non-Ohmic curves skewed 
the opposite way. This difference is due to the response of 
P{u) (the PDF of individual links) to the rise in current. 
In both Ohmic and non-Ohmic regimes the net resistance 
of the system is determined by the largest breaks. But 
in the non-Ohmic case there is almost a hard cutoff « 
ui on the largest possible u (Fig. |9]). In other words, 
breaks with u~> uj are effectively eliminated,^* making 
the large-resistance side of the PDFs of ln(i?/i?o) S> u/ 
drop sharply as well. 

Another result of removing the highly resistant links is 
the PDF's approach to the Gaussian shape. By reducing 
the spread of the link resistances, it brings the system 
closer to the conditions at which the central-limit theo- 
rem is obeyed. This can be seen in Fig. [TO] where the 
curves become narrower and more Gaussian at lower uj. 



Also plotted in Fig. 10 are PDFs obtained by an ap- 
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proximate but much faster method, which utilizes our 
analytical formulas for P{u). We call this the PDF- 
algorithm. The idea is as follows!^ The resistance of 
the system is given by the sum over all links, 

R^Ro^e'^' . (45) 

i=l 

Under the assumption that the link resistances are in- 
dependent random variables, each with the same PDF 
P(u), it can be shown that 

Pu{U) = ^ / exp (C/ - ite^) g{t)dt , (46) 
g{t) = exp / p(w)[exp(ite") - l]du\ , (47) 




\n{RIR^ 



FIG. 10: The PDF of the logarithm of the total resistance 
R. The values of ui are indicated next to each curve. The 
simulation parameters are the same as in Fig. |2] The smooth 
curves are obtained using the PDF algorithm, the markers 
are from the shortest-path simulations. 



This is equivalent to the formulas given by RR in Refs.|3HI 
and llll For convenience of the reader, we include a quick 
derivation. For independent variables the cumulants of 
the sum are equal to the sum of the cumulants To cal- 
culate the latter we notice that the number of breaks of 
size (u, u -|- du) has the average value dN{u) — Lp{u)du. 
The actual number is random and has the Poisson distri- 
bution. Therefore, its contribution to nth cumulant of 
R/Ro is e^'^dN. The total cumulant is 



L / p{u)e''''du. 



(48) 



Reconstructing the characteristic function Q{t) from the 
cumulants in a standard way,'^we obtain Eq. ( [47| . Tak- 
ing its Fourier transform and making the change of vari- 
able from i? to [/, we recover Eq. (46 1. 



Certainly, the resistances of the links are not truly un- 
correlated; however, since R is dominated by the largest 
breaks, which are rare and well-separated, this should be 
a good approximation. Note that in Ref. [37] an attempt 
was made to include correlations between adjacent links. 
As mentioned above, it does not compare well with our 
simulations. 

In practice, even a numerical integration of the 
strongly oscillating functions in Eqs. (46 1 and (47 1 is diffi- 
cult. We found it easier to directly implement Eq. (45 ) in- 



stead. To this end we draw Ui from the distribution P{u) 
using a Monte-Carlo sampling (the usual acceptance- 
rejection algorithm). After [Eq. ( [41] )] of such resis- 
tances are generated, the total resistance of the wire is 
obtained by summing them. Figure [TO] illustrates that 
the PDFs obtained from the shortest-path simulations 
and from the PDF-algorithm are in a good agreement. 
The curves produced by the latter are much more smooth 
because we could apply it to a larger number of disorder 
realizations: 10^. 



V. CONDUCTANCE-VOLTAGE 
CHARACTERISTICS 

Having studied the statistics of individual hops that 
contribute to the ID transport, we can now move to the 
analysis of macroscopic transport properties. In exper- 
iment, such transport properties are measured either as 
a function of current or as a function of voltage. In the 
former case, the ensemble averaging gives the average re- 
sistance (i?); in the latter — the average conductance 
(G). If a large number of nominally identical wires is 
available simultaneously, this can be done in a single mea- 
surement, connecting them, respectively, in series and in 
paralleLiJ^ Otherwise, one can try to create the members 
of an ensemble one by one by varying gate voltage or 
other parameters of a single wire.*^ 

Since our shortest-path algorithm is formulated at a 
constant current (i.e., constant uj), one may naively 
think that it is able to provide only the distribution of 
resistances. This is not so. Let us show that the PDFs of 
conductances and resistances are uniquely related even 
in the non-Ohmic regime. 

Denote the PDF of having a given total voltage F at a 
fixed current / by Pv{V\I) and the PDF of having a given 
current / at a fixed total V by Pi{I\V). By inspecting 
the V-I curves sketched in Fig. |11[ we can write down 
the following continuity equation: 

^jPv{V\I) + ^Pi{I\V)^i). (49) 
Integrating with respect to voltage, we get 

V 

P,{I\V) = -^^ J Pv{V'\I)dV'. (50) 



As an application, let us show how the average conduc- 
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I+dl 



FIG. 11: A sketch of V~I curves for an array of different wires. 
The continuity equation ( 49 1 follows from the conservation of 



the number of curves piercing the differential area element 
bounded by the dashed lines. 



tance Gy at a given fixed voltage V , 



(51) 



can be calculated. 



In view of Eq. ( 50 1 , Gy can also be written as 



oo V 

Gv = - T-f^J Pv{V'\I)dV'. (52) 











We integrate this by parts and change the notation for 
the measure in the second integral from Pv{V'\I)dV' to 
Pii{R'\ui)dR' . We arrive at the formula 

-oo 

for the desired average conductance at a fixed voltage. It 
is easy to see that in the Ohmic limit, V ^ Eq. (53) 



coincides with the average conductance at a fixed current, 
J PR{R'\oo)dR'/R', as expected. 

To evaluate Gy as a function of V one needs to know 
Pii{R'\uj). We obtained it by the following procedure. 
We divided the interval of V we are interested in into 
a number of bins. We took an interval of ui from 5 to 
about ui — 3um and in turn divided it into equidistant 
steps ui{j), I < j < Nj = 1000, spaced by Auj. For 
each ui{j) we generated Ngam = 200 samples, i.e., sets 
of individual u's, drawn from the distribution P{u) 
using the acceptance-rejection algorithm. We converted 
the integrals in Eq. (53 1 into discrete sums. 



G\ 
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Rq 



(54) 

where Ri{j) is the total resistance of ith set for a given 
j, and then evaluated them numerically. 

The simulations were done for um = 5, 7.5, 10, 12.5, 
and 15. The control parameter was T while all other 



values — a, g, L, and Nu — remained the same. Later 
we realized that in the non-Ohmic regime the number of 
hops gradually increased with current. Equation (41 1 
remains accurate only for uj > um- Therefore, only 
uj > Um points were included when plotting the five 
curves in Fig. [3) 

Alternatively, Gy can be reduced to a numerical 
quadrature, which this time contains no oscillating in- 
tegrands. This is possible because Gy is dominated by 
large conductances, for which the saddle-point approxi- 
mation in Eq. (46 1 is legitimate. After a straightforward 
derivation, one obtains 



Gv _T 1 dui J^fVe^' ^ 

-oo 




X ^ / — exp ( Jit + Jt))dt , 
27r 



(55) 



Jn^N^J P{u) [exp{nu - te") - S^^o] du , (56) 



where n = 0, 1,2, and Sij is the Kronecker symbol. All 
these integrals are rapidly converging, so that their nu- 
merical evaluation should cause no difficulty. However, 
we deemed the quality of the curves shown in Fig. [3] suffi- 
cient [these curves were obtained from Eq. (54|]. There- 
fore, we did not pursue this alternative method. 



VI. DISCUSSION 

At this point, let us recapitulate our findings. To the 
best of our knowledge we presented the first reliable cal- 
culation of the statistics of resistances in ID VRH net- 
work, both in Ohmic and non-Ohmic regimes. Compar- 
ing with the previous theoretical work, we showed the 
importance of the correction to the PDF P{u) proposed 
in Ref. 36 We demonstrated that without this "sublead- 
ing" term the conductance could be significantly overes- 
timated, see Fig. [T] Figure |6] further illustrates the im- 
portance of such corrections by showing that there are no 
obvious diamond-like or hexagonal voids in the energy- 
position space invok ed in the derivations of the leading 
asymptotic behavior .ISIill 

Next, our calculations have verified the earlier analyti- 
cal predictions^^ that large breaks are progressively elim- 
inated at higher voltage, and that the PDF of resistances 
becomes more narrow, see Fig. |10[ This disappearance 
of highly resistive hops equalizes different samples, mak- 
ing the averages of parallel and series setups of the wires 
approach the same value. 

Let us now turn to experiments. Unfortunately, we 
could not find a clear evidence of the predicted behav- 
ior in published literature. A dedicated experiment to 
probe mesoscopic conductance fluctuations in non-Ohmic 
regime is desired as it was not on the agenda in previous 
studies of ID VRH. At least two other caveats must also 
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be kept in mind. First, most of "ID" electron systems 
studied experimentally were not truly one-dimensional. 
They either consis ted of many parallel chain^^l^ or had 
multiple subband^^^m or were bulk samples with a 
large aspect ratio.'^^'^Uguch systems may behave as effec- 
tively ID but only at low enough T . Finally, our model 
of disorder where LS are treated as points in the energy- 
position space may or may not be relevant for some of 
these experiments (see more below). 

Turning to some specific examples, we consider first 
the measurements done on polydiacetylene single crys- 
tals,'32l which are quasi-lD materials. The Ohmic trans- 
port is consistent with ID VRH behavior, showing a 
crossover from a simple exponential at relatively high 
temperatures, InG w —Ah/2T, to a stretched expo- 
nential InG « -(A;/r)T with 7 = 0.5-0.75 at low T. 
Just as in our simulations, there is a substantial differ- 
ence between A/j and A;. For instance, in sample SI 
Ah = 320 K and A/ = 2570 K w 8A,,. In the same 
sample at high electric fields Eq. Q is observed, with 
8To/a = 0.049 eV/nm (in our notations). Assuming that 
Tq « Ah, this gives a reasonable estimate of the local- 
ization length a — 4.3 nm. At modest fields, the trans- 
port data were fitted to Eq. ^ and Lc was extracted. 
It was seen to have the same temperature dependence 
Lc oc T""'^, as in our simulations. Moreover, the numer- 
ical value of Lc is close to what we find. For example, 
Lc = 32.5 nm at T = 25 K in t he expe riments,!^ which 
can be compared to Lc ^ 1.9a-\/2To/r = 40 nm that we 
find, cf. Fig. |4] 

Next, let us consider another experiment, which was 
done on arrays of GaAs quantum wires. The depen- 
dence of G on F and T that we have calculated here is in 
a reasonable agreement with some of those experimental 
results but some strong deviations are also apparent. For 
example, in the simulations the range of activated behav- 
ior in the Ohmic regime spans at best two decades in G. 
In the experiment, it is much wider (three decades), and 
occupies most of the temperature range T > 0.2 K where 
G was reported. We were able to fit the experimental 
G(0,T) only by imposing a rather strong power-law de- 
pendence of the prefactor: Go c>c T^-^. From such a fit 
we obtained Tq = 6.2 K (in our notations). 

In the non-Ohmic regime, the initial rise of G with F 
is again exponential over approximately one decade, see 
Fig. [12] However the behavior of parameter Lc in this 
exponential law was deemed to be surprising in Ref. 1171 
Therefore, let us discuss it. Physically, Lc is the distance 
between "critical hops" in a sample, i.e., those highly 
resistive links that generate the dominant portion of the 
total voltage. In a typical sample, length Lc has to be 
much larger than the average hop length uj\/a. In fact, at 
low T one would naively expect Lc to be of the order of 
the sample length L. This is because in a typical sample 
all the voltage drops on a single break. At higher T, 
where the activated trans por t is observed, the voltage 
is shared by many breaks,'^ and so Lc is supposed to 
decrease exponentially. However, this is not what was 
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FIG. 12: Conductance of sample 1 of Ref. [T7| as a function 
of the scaled electric field Fa/To (markers). Here To — 6.2 K 
is determined from the best fit of Eq. ^ to the Ohmic con- 
ductance (not shown) and a = 0.4 pm. Temperature (in K) 
is indicated next to each data set. The best fits to Eq. Q 
are shown by the lines. The rightmost curve is Eq. Q. The 
prefactor Go is chosen such that the relation between Eq. Q 
and the uppermost data trace (corresponding to um ~ 5) is 
similar to that in Fig. [S] 



observed. At low T, two out of three samples measured 
in Ref. 17 had Lc ~ L/50, while Lc of the third was 
about L/10. As T was increasing, Lc was decreasing but 
rather slowly, perhaps, as T^^/^. 

In light of our findings, this behavior of Lc is not sur- 
prising. The above reasoning does not take into account 
that the measurements were done not on a single wire but 
on several hundreds of them, connected in parallel. It is 
logical to assume that some wires conducted much bet- 
ter than others because they happened to have no breaks. 
These wires could short out the wires which were poor 
conductors, reducing the net Lc down to the typical hop- 
ping length. 

We now demonstrate explicitly that Lc extracted from 
our model has numerical values and functional behavior 
similar to what was measured experimentally. Our Lc, 
which was found by fitting the low-voltage part of G{F) 
curves m Fig. [3] to Eq. ^ is plotted in Fig. g The 
intervals of Tq [T are different in our simulation and the 
experiment; however, there is a small overlap. For our 
leftmost point, TqIT — 12.5 we have LjLc ~ 30, similar 
to the numbers quoted above. 

The problem arises when we consider the high-field be- 
havior reported in Ref. '17', Experimental G(F, T) curves 
tend to approach a common T-independent limit, as in 
our calculations. Fig. [3] However, this limit is strongly 



underestimated by our Eq. (|4|), see Fig. [12j While we 
do not know the origin of this discrepancy, it is possible 
that the different behavior seen in the two experiments is 
just another example of a dilemma, which has a long his- 
tory in the VRH literature. Previously, it was discussed 
mostly in the context of bulk materials where majority of 
experiments have been done so far. However, it is tempt- 
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ing to make a comparison with our ID case because the 
VRH exponent of the Efros-Shklovskii law in any dimen- 
sion nominally coincides with the ID Mott law exponent 
7 = 1/2. 

The essence of the dilemma is as follows. There are a 
number of systems where non-Ohmic behavior does fol- 
low Eqs. ^ and Q that we have observed in our simula- 
tions. However, this is usually the case when parameter 
A in Eq. ([!]) is large, say, tens or hundreds of K. Very dif- 
ferent and still poorly understood behavior occurs when 
Tq is relatively small (according to one study,!^ when 
\/TqIT < 12). The high-field nonlinearities in this sec- 
ond group are much stronger. In the extreme cases , the 
I-V characteristic was determined to be S'-shaped,'^^'^ 
which led to hy steretic conductivity j umps by orders of 
magnitudd^^l^S and circuit oscillations.'^^Hfil Interestingly, 
in systems that show conductivity jumps the Ohmic con- 
ducta nce sh ows a simple activation rather than VRH be- 
havior.l^SES! 

It has become commorJ^U^' 45 | 47 | 48 | 49 | 50 | 5i | 5 2] ^ 

tribute strong nonlinearity and S'-shaped I-V to elec- 
tron overheating. It is assumed that G is the function of 
the electron temperature Tg, which can be much higher 
than the ambient temperature T. A phenomenological 
equation is postulated. 



Q ^ GF^ ^ a{T^ - T^) , 



(57) 



where a and (3 are adjustable constants. (Usually, 4 < 
(3 < 8.) This equation is supposed to represent the 
balance between the Joule heat delivered into electron 
system from the external field and the heat transferred 
from electrons to phonons. Surprisingly, this equation 
has been shown to provide an accurate description of 
some VRH systems, including the the one we are trying 
to make comparison to.'^^'^ 

By itself, the idea of hot electrons is not objection- 
able. Actually, our Eq. ^ can be viewed as the ID 
Mott law with the electron temperature Tg ~ Fa (sim- 
ilar to Refs. 1^ and "531. The difficulty is that the re- 
quired Te is unusually large. Indeed, let us define the 
length Lc-ph = T^/F. It has the physical meaning of a 
characteristic distance over which an electron must be 
accelerated by the external field to gain the extra energy 
Te 3> T. In our model, where LS are treated as points, 
the largest achievable ie-ph is of the order of a. Electrons 



cannot propagate farther without suffering an exponen- 
tial decay. Yet to get a stronger I-V nonlinearity than 
predicted by our Eq. Q, ic-ph must exceed a. For exam- 
ple, to reproduce the high-field part of the data shown in 
Fig. 12 we need perhaps io-ph ^ 10a. 

In principle, Lc-ph 3> a is possible if the disordered 
system is a granular metal or equivalently, an array of 
random- sized quantum dots. In this case the upper 
bound on Le-ph is presumably set by the size of metal- 
lic grains, while the exponential decay length a is much 
smaller, being suppressed by weak tunneling between the 
grains. The granular-metal model can also explain a wide 
range of the activated Ohmic behavior as a manifestation 
of the Coulomb blockade. Finally, it has been suggestecP^ 
that the conductivity jumps may be related to lifting of 
the Coulomb blockade by collective depinning. Trans- 
port in a ID version of this model was recently studied 
in a paper co-authored by one of u^^ but the case of 
extremely strong fields was not considered. It remains to 
be seen whether this model can yield a better agreement 
with the experiments. 

It has been frequently speculated that the overheating 
is driven by the electron interactions, which we did not 
address here. The simplest way to introduce some inter- 
action effects into the existing formalism is to consider 
larger dielectric constant e > 1. The importance of such 
effects requires further study. 

Finally, as mentioned above, most of electron sys- 
tems studied should behave as effectively ID only at low 
enough T. The dimensional crossover as a function of 
temperature in a strip geometry has been studied by RR 
in Ref. 55, It would be interesting to investigate the 
electric-field counterpart of this crossover. 

In conclusion, we showed that numerical simulations 
such as those we carry out in this paper can serve as a 
valuable tool in studying VRH transport. We hope that 
our results would stimulate further experimental work 
on both "conventional" (semiconductor wires) and novel 
(nanotubes, nanofibers, graphene ribbons) ID and quasi- 
ID materials. 
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